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Abstract 

Many problems in physics, material sciences, chemistry and biology can 
be abstractly formulated as a system that navigates over a complex energy 
landscape of high or infinite dimensions. Well-known examples include phase 
transitions of condensed matter, conformational changes of biopolymers, and 
chemical reactions. The energy landscape typically exhibits multiscale fea- 
tures, giving rise to the multiscale nature of the dynamics. This is one of the 
main challenges that we face in computational science. In this report, we will 
review the recent work done by scientists from several disciplines on probing 
such energy landscapes. Of particular interest is the analysis and computa- 
tion of transition pathways and transition rates between metastable states. 
We will then present the string method that has proven to be very effective 
for some truly complex systems in material science and chemistry. 
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1. Introduction 

Many problems in biology, chemistry and material science can be formulated 
as the study of the energy or free energy landscape of the underlying system. Well- 
known examples of such problems include the conformational changes of macro- 
molecules, chemical reactions and nucleation in condensed systems. Very often the 
dimension of the state space is very large, and the energy landscape exhibits a hi- 
erarchy of structures and scales. These problems are becoming a major challenge 
in their respective scientific disciplines and are beginning to receive attention from 
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the mathematics community. In this article, we report recent work in this direction. 
For a detailed account, we refer to [pj |6|, [zj, 

We begin with a simple example. Plotted in Figure 1 is the solution of the 
stochastic differential equation 



dx{t) = -V^V{x{t))dt + y/sdW{t) 



(1.1) 



where the potential 



and dW{t) is Gaussian white noise, e 
perturbation, the solution would be x{t) 



■ 0.06, a;(0) 
x{0) = -1. 



(1.2) 

- —1. Without the random 
Indeed the deterministic part 



of dynamics in (1.1) does nothing but taking the system to local equilibrium states. 



With the random perturbation, the solution, over long time, exhibits completely 
different behavior. It fluctuates around the two local minima oi V,x = —1 and 1, 
with sudden transitions between these two states. The time scale of the transition, 
tM is much larger than the time scale of the fluctuation around the local minima, 
tR. For this reason, we refer to a; = —1 and 1 as the metastable states. 



eps = 0.06 






Figure 1. Time series of the solution to the stochastic differential 
equation (1.1), with e = 0.06. 



Obviously the transition between the metastable states is of more interest 
than the local fluctuation around them. The transition time is much larger since it 
requires the system to overcome the energy barrier between the two states. This is 
only possible because of the noise. When e is small, a huge noise term is required 
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to accomplish this. For this reason, such events are very rare, and this is the origin 
of the disparity between the time scales i*/ and tji. 

This simple example illustrates one of the major difficulties in modeling such 
systems, namely the disparity of the time scales. It does not, however, illustrate 
the other major difficulty, namely, the large dimension of the state space and the 
complexity of the energy landscapes. Indeed for typical systems of interests the 
energy landscape can be very complex. There can be a huge number of local minima 
in the state space. The usual concept of hopping over barriers via saddle points 
may not apply (see |^). 

In applications, the noise comes typically from thermal noise. In this case, 
we should note that even though the potential energy landscapes might be rough 
and contain small scale features, the system itself experiences a much smoother 
landscape, the free energy landscape, since some of the small scale features on the 
potential energy landscape are smoothed out by the thermal noise. 

Our objective in modeling such systems are the following: 

1. Find the transition mechanism between the metastable states. 

2. Find the transition rates. 

3. Reduce the original dynamics to the dynamics of a Markov chain on the 
metastable states. 

Our discussion will be centered around the following model problems: 

7i(i) = -VVix{t)) + ^W{t) (1.3) 



or 

mx{t) + 7i(i) = -VV{x{t)) + y/^W{t) (1.4) 
e is related to the temperature of the system by e = where fcs is the 



Boltzmann constant. We refer to ( |l.3[ ) as type-I gradient flow and ( |l.4D as type-II 
gradient flow. 

Before proceeding further, let us remark that there is a very well-developed 
theory, the large-deviation theory, or the Wentzell-Freidlin theory that deals 
precisely with questions of the type that we discussed above. However as was 
explained in ||^, this theory is not best suited for numerical purpose. Therefore 
we will seek an alternative theoretical framework that is more useful for numerical 
computations. 



2. Transition state theory 

Transition state theory (TST) has been the classical framework for ad- 
dressing the questions we are interested in. It assumes the existence and explicit 
knowledge of a reaction coordinate, denoted by q, that connects the two metastable 
states. In addition it assumes that along the reaction coordinate there exists a well- 
defined transition state, which is typically the saddle point configuration, say at 
q = 0, and the two regions {q < 0} and {q > 0} defines the two metastable regions 
A and B. For these reasons, transition state theory is restricted to cases when the 
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system is simple and the energy landscape is smooth, i.e. the energy barriers are 
larger than the thermal energy ksT. 

Knowing the transition state, TST calculates the transition rates by placing 
particles at the transition state, and measuring the flux that goes into the two 
regions. For example, the transition rate from A to i? is given approximately by 

kA^B = ^ I mSr{qitmq{t))dfiAiqm (2.1) 



Zo 
where 

Zo = JdfiAiqm (2.2) 

Here Sr is the surface delta function at g = 0, 6' is the Heaviside function, ^a is the 
Lebesgue measure restricted to A. For a system with a single particle of mass m 
and potential V, this gives 

kA^B = -^e-^ (2.3) 
where SE is the energy barrier at the transition state, ujq = m^"* ) ' denotes 



the location of the local minimum inside A. Formulas such as (2.3) are the origin 
of the Arrhenius law for chemical reaction rates and Boltzmann factor for hopping 
rates in kinetic Monte Carlo models. 



3. Reduction to Markov chains on graphs 



For simplicity, we will discuss mainly type-I gradient flows ( |l.3[) . The Fokker- 
Planck equation can be expressed as 

|,,„.V.(„„V(^)) ,3.) 

where ps is the equilibrium distribution 

1 

Ps{x) = — e ^B-r 

Z is the normalization constant Z = J^„ e ''b^ dx. The states of the Markov chain 
consist of the sets {Bj}j^i, where satisfies 

1. The i?j's are mutually disjoint 
2. 



Psix)dx = I + o{kBT) (3.2) 

where B — Uj^iBj. An illustration of the collection the is given in Figure 

2. {Bj} depends on T. As T decreases, the B/s exhibits a hierarchical structure. 
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Having defined the states of the Markov chain, we next compute the transition 
rates between neighboring states. Denote by A and B two such neighboring states. 
We would hke to compute the transition rate from A to B. Without loss of gener- 
ality, we may assume J = 2. Let B2 be the metastable region containing A and 
B respectively, and let nj{t) = p{x,t)dx, Nj = ps{x)dx, j = 1,2. Applying 
Laplace's method to (3.1) we get ^ 

. , , £ f n2{t) ni{t)\ 

nj[t) — — I — — — — 1 + higher order terms (3.3) 



where 



N2 Ni 



SO{a) 





^"(a)| (3.4) 



{(p'~'{a), < a < 1} is a so-called minimal energy path, to be defined below, {S'~'{a)} 
is the family of hyperplanes normal to (p^. 

The minimal energy path (MEP) is defined as follows. If V is smooth, then 
ip° is a MEP if 

(Vy)^(^°(a))=0 (3.5) 

for all a G [0, 1], i.e. WV restricted to is parallel to ip'^. In general there is not a 
unique ip^ that is particularly significant, but rather a collection (a tube or several 
tubes) of paths contribute to the transition rates. However, one can define a MEP 
self-consistently via the equation 

(^O(a) = [ xe^^dx (3.6) 

Z[aj J so (a) 
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where Z{a) — Jgo^^^ e ''b'^ dx. In the case when V has two scales: V — V + 5V , 
\5V\ < OiksT), and V is smooth, then ^3° can be defined as the MEP of V. 

For type-I gradient systems, if the two metastable sets are separated by a 
single saddle point, then the MEP is the unstable manifold associated with the 
saddle point. MEP for type-II gradient systems is less trivial. In this case (3.3)- 
(3.6) have to be modified ||]. Consider the simple example 



q =p 



where V{q) = j{l — q"^)^. It has two local equifibrium states A = (—1,0) and 
B = (1,0). The MEP that connects A and B is plotted in Figure 3. It is not a 
smooth curve. The velocity is reversed at the saddle point. To verify that the MEP 
does reflect the true behavior of the transition path, we also plot the transition path 
obtained from direct simulation of the stochastic differential equation. 




Figure 3. MEP for type-II gradient systems. The red line is 
the MEP in phase space, the green line is the transition path 
computed from solving the stochastic differential equation. 



MEP is a very important concept since it defines the "most probable" transi- 
tion path from which transition rates can be computed via equation (3.3). However 
we should emphasize that from a numerical point of view our task is not that 
of a conventional optimization or control problem, since there is not an objective 
function that we can easily work with. Instead our aim is to perform importance 
sampling in path space to sample the paths that contribute significantly to the 
switching. 
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Finally, if there exists a MEP that connects two metastable sets without going 
through a third one, then we connect these two metastable sets by a link. In this 
way, we form a graph. The original dynamics is then reduced to a Markov chain on 
this graph. 

4. Previous numerical techniques 

A variety of numerical techniques have been developed, most prominently 
in chemistry, but also in biology and material science, for computing MEPs and 
sometimes transition rates. Among the most well-known techniques in the chemistry 
literature are the nudged elastic band method and the transition path sampling 
technique. The nudged elastic band method (NEB) aims at computing the 
MEP defined by (3.5). It represents the MEP by a discrete chain of states. These 
states evolve according to the potential forces of the system. To prevent the states 
from all falling to the two local equilibrium states, a spring force is applied to 
neighboring states to penalize the non-uniformity in the distribution of the states 
along the chain. This by itself may cause convergence to a path which is not a 
MEP. Hence a nudging technique is used, namely only the normal component of 
the potential force and the tangential component of the spring force is applied. 
NEB is a very effective method for small systems with relatively smooth energy 
landscapes. It has two main drawbacks. One is that it is highly inefficient and may 
not even be applicable to systems with rough energy landscapes. The other is the 
choice of the elastic constant. A large elastic constant requires a small time step 
in the evolution of the states. A small elastic constant will not achieve the desired 
uniformity of the states and hence will not give the required accuracy for the energy 
barrier. 

A second important technique is transition path sampling (TPS) H, js). This 
method aims at complex systems with rough energy landscapes by developing a 
Monte Carlo technique that samples the path space. Its efficiency hinges on the 
ability to produce new accepted paths from old ones. 

Other techniques include the ridge method, blue mooth sampling, etc. 
Often these methods require knowing beforehand the reaction coordinate. 

Elber et. al propose to minimize the Onsager-Machlup action as a way of 
finding the most probable path for macromolecular systems jl^. The Onsager- 
Machlup action is the same as the Wentzell-Freidlin action. From a numerical 
point of view, there are certain difficulties associated with minimizing this action 
functional. These issues are discussed in Q. 

5. The string method 

The basic idea in the string method, developed in [§, ^, 0], is to represent 
transition paths by their intrinsic parameterization in order to efficiently evolve and 
sample paths in path space. It has two versions. The zero temperature version is 
designed for smooth energy landscapes. The finite temperature version is designed 
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for rough energy landscapes in which case thermal noise acts to smooth out the 
small scale features. 

The simplest example of a zero-temperature string method is to evolve curves 
in path space by the gradient flow 

- -i\7V)^iip{a,t)) +ria,t)f{a,t) (5.1) 

Here f is the tangent vector of the curve {tp{-,t)}, {WV)'^{(p) denotes the compo- 
nent of W normal to f, r is the Lagrange multiplier that enforces certain specific 
parameterization of the curves. For example if we require equal arclength parame- 
terization, then we need ^ItySa] =0, i.e. 

/• 1 /•a 

r{a,t)=a \/V{ip{a' ,t)) ■ T{a' ,t)da' ~ \/V{ip{a' ,t)) ■ f{a' ,t)da' (5.2) 



We call such curves with intrinsic parameterization strings. 

In practice the strings are discretized into a collection of points. These points 
move according to the normal component of the potential force. A reparameteriza- 
tion step is applied once in a while to enforce the proper parameterization of the 
strings. 

The finite temperature string method is designed for systems with rough energy 
landscapes, particularly the case when the potential can be expressed in the form 

V{x) = V{x)+dV{x) 

where V is smooth and \6V\ < 0[kBT). In this case we would like to compute the 
MEP of V without first computing V explicitly. This is achieved by creating an 
ensemble of a special type. Our computational object will be a string connecting 
the two metastable sets, together with a family of probability measures on the 
hyperplanes normal to the string. Consider the stochastic equation 

^r(a,t) - -P:?:(Vl/(^"(a,t))) +rO(a,t)f"(a,t) +P^;-(a,i) (5.3) 

where rj^ is Gaussian noise with mean and correlation 



E,7"(a,t),7"(a',r) 



2kBTS{t~T), ifa = a' 
0, if a 7^ a' 



The projection operator Pj^ is defined by projecting to the hyperplane normal to 
the string {(^°(a,i),0 < a < 1}, where 

^°(a,i)=E(p"(«,t) 

f*^(Q!, t) is the tangent vector of ip^ at a, r° is the Lagrange multiplier that enforces 
proper parameterization of ip'^ . 
Theorem 1.1. 
1. The statistical steady state of (5_^) satisfies: 

n If _ZM 

V \.ot) = ——l xe "B^dx (5.4) 
^(") Js°(a) 

where S'^{a) is the hyperplane normal to ifP at a, Z(a) — /50(q,) g ''b'^ dx. 
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2. The stationary distribution of (5.3) is given by the family of distributions 



Z{a) 



'5s"(a){x) 



(5.5) 



Knowing {(/s"} and {/ia}, the transition rates and free energy landscapes can 
be computed, see 0. 

The finite temperature string method is applied to the perturbed Mueller 
potential 

V{x)^V^{x)+5V{x) 

where Vm{x) is the so-called Mueller potential (see 5V is a random pertur- 

bation. The results of the (finite temperature) string method is shown in Figure 4. 
Also plotted is the MEP of Vm (since V = Vm is explicitly known for this particular 
example) as well as the fluctuations around (/3°. 



1.5 



0.5 




Figure 4. Effective MEP and local fluctuations for the perturbed 
mueller potential. The red curve is the MEP for V — Vm- The 
black curve is the MEP computed from the finite temperature 
string method. The green curves show the size of the fluctuations. 



6. Concluding remarks 

There are several important topics that we did not cover in this brief report. 
These include the effect of dynamics, non-gradient systems, and acceleration tech- 
niques. These are discussed in |Q, |^, ^ |T^. Also found in these references are 
applications of the ideas discussed here to thermal activated reversal of magnetic 
thin films, models of martensitic transformations, and the formation of Cgo from 
60 carbon atoms. The last example is a case when the barrier is entropic. Even 
though the potential energy is mainly going downhill, the free energy has barriers 
because of entropic effects. Such examples are found frequently in biopolymers. 
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Prom a numerical point of view, our main idea for overcoming the difficulty 

caused by the disparity of the times scales is to reformulate the problem as a bound- 
ary value problem instead of initial value problem, since we have some knowledge of 
the initial and final state of the system. Compared with other existing methods that 
assume explicit knowledge of a reaction coordinate, our method finds the reaction 
coordinate self-consistently during the computation. 

The topic discussed here is relatively new in applied mathematics, but it is of 
paramount importance in science and particularly computational science. Progress 
in this area will likely have a fundamental impact in many areas of applications. 
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